!! Immersed Boundary module
module immersed_boundary
use omplib
use technical
use general
use mod_2dflu
implicit none
save
!
integer :: nboundary_particle=1024
integer :: nneighbour=0
double precision,allocatable,dimension(:)::Xboundary,Yboundary,dell,Vxboundary,Vyboundary
double precision,allocatable,dimension(:)::Xboundary_fixed,Yboundary_fixed
double precision,allocatable,dimension(:,:)::fboundary
double precision,allocatable,dimension(:,:,:)::FEuler
integer,allocatable,dimension(:,:) :: neighbour_coord_list
double precision :: elastic_constant=10.
logical :: lmoving_boundary=.false.
character (len=labellen) :: immersed_boundary_type='none'
double precision :: snake_frac_width=0.5,snake_amp=0.1
integer :: snake_wavenumber=2, dirac_cutoff_grid=2
double precision :: dirac_delta_cutoff 
double precision :: dq
!
private
public :: read_immersed_boundary,allocate_immersed_boundary,initialize_immersed_boundary,evolve_immersed_boundary
public :: add_boundary_force,get_immersed_boundary_rhs
namelist /immersed_boundary_pars/ &
  nboundary_particle,elastic_constant,immersed_boundary_type,&
  snake_frac_width,snake_amp,snake_wavenumber,dirac_cutoff_grid
contains
!**********************************************************!
subroutine read_immersed_boundary(unit,ierr)  
  integer, intent(in) :: unit
  integer, intent(inout), optional :: ierr
!----------
  if (present(ierr)) then
    read(unit,NML=immersed_boundary_pars,IOSTAT=ierr)
  else
    read(unit,NML=immersed_boundary_pars)
  endif
  rewind(unit)
endsubroutine read_immersed_boundary
!!--------------------------------------
subroutine allocate_immersed_boundary()
  allocate(Xboundary(nboundary_particle),Yboundary(nboundary_particle))
  allocate(dell(nboundary_particle))
  allocate(Vxboundary(nboundary_particle),Vyboundary(nboundary_particle))
  allocate(Xboundary_fixed(nboundary_particle),Yboundary_fixed(nboundary_particle))
  allocate(fboundary(nboundary_particle,2))
  allocate(FEuler(n1,n1,2))
endsubroutine allocate_immersed_boundary
!!--------------------------------------
subroutine initialize_immersed_boundary()
  dirac_delta_cutoff=dx*dfloat(dirac_cutoff_grid)
  write(*,*) ' Using dirac_cutoff_grid=',dirac_cutoff_grid
  write(*,*) ' Using dirac_delta_cutoff=',dirac_delta_cutoff
  call make_neighbour_coordinates_list()
  select case(immersed_boundary_type)
    case('snake')
      call set_immersed_boundary_snake()
    case('peristalsis')
      call set_immersed_boundary_peristalsis()
    case default 
      call fatal_error('initialize_immersed_boundary:','please select immersed_boundary_type')
    endselect
  Xboundary_fixed=Xboundary
  Yboundary_fixed=Yboundary
!
endsubroutine initialize_immersed_boundary
!!--------------------------------------
  subroutine make_neighbour_coordinates_list()
    integer :: ix,iy
    double precision :: rx,ry,rr
    integer :: nnumber,ineighbour
      nnumber=0
      do ix=-dirac_cutoff_grid-1,dirac_cutoff_grid+1;do iy=-dirac_cutoff_grid-1,dirac_cutoff_grid+1
        rx=dfloat(ix)*dx; ry=dfloat(iy)*dy
        rr=sqrt(rx*rx+ry*ry)
        if (rr .le. (dirac_delta_cutoff+dx) ) then
          nnumber=nnumber+1
        endif 
      enddo;enddo
      write(*,*) 'Each grid point has  ',nnumber,'neighbours, allocating accordingly'
      allocate(neighbour_coord_list(nnumber,2))
      nneighbour=nnumber
      nnumber=0
      do ix=-dirac_cutoff_grid-1,dirac_cutoff_grid+1;do iy=-dirac_cutoff_grid-1,dirac_cutoff_grid+1
        rx=dfloat(ix)*dx; ry=dfloat(iy)*dy
        rr=sqrt(rx*rx+ry*ry)
        if (rr .le. (dirac_delta_cutoff+dx) ) then
          nnumber=nnumber+1
          neighbour_coord_list(nnumber,1)=ix 
          neighbour_coord_list(nnumber,2)=iy 
        endif 
      enddo;enddo
      write(*,*) 'Now writing the list to file ..'
      open(unit=5,file='neighbour.list',status='unknown')
      do ineighbour=1,nnumber
         write(5,*)neighbour_coord_list(ineighbour,1), neighbour_coord_list(ineighbour,2)
      enddo
      close(5)
  endsubroutine make_neighbour_coordinates_list
!!--------------------------------------
subroutine set_immersed_boundary_snake()
  integer :: ibp
  double precision :: dbp,snake_kk,snake_width,xbp,ymid,dXb,dYb
  lmoving_boundary=.false.
  snake_width=snake_frac_width*length
  dbp=length/(nboundary_particle/2)
  snake_kk=dfloat(snake_wavenumber)
  ymid=(ymax+yy0)/2.
  do ibp=1,nboundary_particle/2
! upper boundary of the snake
    xbp=xx0+(ibp-1)*dbp
    Xboundary(ibp) = xbp
    Yboundary(ibp) = ymid+snake_width/2.+snake_amp*sin(snake_kk*xbp)
    dXb=1.
    dYb=snake_amp*snake_kk*cos(snake_kk*xbp)
    dell(ibp) = dbp*sqrt(dXb*dXb+dYb*dYb)
    dell(nboundary_particle/2+ibp) = dell(ibp);
! lower boundary of the snake
    Xboundary(nboundary_particle/2+ibp) =  xbp
    Yboundary(nboundary_particle/2+ibp) = ymid-snake_width/2.+snake_amp*sin(snake_kk*xbp)
  enddo
  call write_input_boundary_points()
endsubroutine set_immersed_boundary_snake
!!--------------------------------------
subroutine set_immersed_boundary_peristalsis()
  integer :: ibp
  double precision :: dbp,snake_kk,snake_width,xbp,ymid,dXb,dYb
  write(*,*) 'Setting up immersed boundary for peristalsis'
  lmoving_boundary=.true.
  snake_width=snake_frac_width*length
  dbp=length/(nboundary_particle/2)
  snake_kk=dfloat(snake_wavenumber)
  ymid=(ymax+yy0)/2.
  do ibp=1,nboundary_particle/2
! upper boundary of the snake
    xbp=xx0+(ibp-1)*dbp
    Xboundary(ibp) = xbp
    Yboundary(ibp) = ymid+snake_width/2.+snake_amp*sin(snake_kk*xbp)
    dXb=1.
    dYb=snake_amp*snake_kk*cos(snake_kk*xbp)
    dell(ibp) = dbp*sqrt(dXb*dXb+dYb*dYb)
    dell(nboundary_particle/2+ibp) = dell(ibp);
! lower boundary of the snake
    Xboundary(nboundary_particle/2+ibp) =  xbp
    Yboundary(nboundary_particle/2+ibp) = ymid-snake_width/2.-snake_amp*sin(snake_kk*xbp)
  enddo
  call write_input_boundary_points()
endsubroutine set_immersed_boundary_peristalsis
!!--------------------------------------
subroutine write_input_boundary_points()
  integer :: ibp
  open(unit=11,file='immersed_boundary_points.ascii')
  do ibp=1,nboundary_particle
    write(11,*) Xboundary(ibp),Yboundary(ibp)
  enddo
  close(11)
endsubroutine write_input_boundary_points
!!--------------------------------------
subroutine get_immersed_boundary_rhs(ux,uy,n1,n2)
  integer, intent(in) :: n1,n2
  double precision, dimension(n1,n2), intent(in) :: ux,uy
  integer :: ibp
  double precision :: xxb,yyb,ddx,ddy
  integer :: xpgrid,ypgrid,xgrid_next,ygrid_next,igrid,jgrid
!
! also set the eulerian boundary force array to zero. 
!
  FEuler = 0.
!
! Calculate the RHS of the equation of motion of boundary points
!
!!
!$OMP PARALLEL DO PRIVATE(ibp)  SHARED(ux,uy,Xboundary,Yboundary,Vxboundary,Vyboundary,FEuler)
  do ibp=1,nboundary_particle
    call get_boundary_velocity(ibp,ux,uy,n1,n2)
!    if(ibp.eq.20) then
!       write(*,*) 'DM,ibp,Vboundary',ibp,Vxboundary(ibp),Vyboundary(ibp)
!       write(*,*) 'DM,ibp,Xboundary',ibp,Xboundary(ibp),Yboundary(ibp)
!       write(*,*) 'DM,ibp,Xboundary_fix',ibp,Xboundary_fixed(ibp),Yboundary_fixed(ibp)
!       xxb=Xboundary(ibp); yyb=Yboundary(ibp)
!       call get_grid_loc(xxb,n1,dx,xpgrid,xgrid_next,ddx)
!       call get_grid_loc(yyb,n2,dy,ypgrid,ygrid_next,ddy)
!    endif
!
! Now calculate the force on the fluid due to the particles. 
! set the force in Eulerian grid to zero to begin with
!
    call calc_boundary_force_lag(ibp)
!    if(ibp.eq.20) then
!       write(*,*) 'DM,ibp,fboundary_lag',ibp,fboundary(ibp,:)
!    endif
    call from_lag2euler(ibp,fboundary(ibp,:))
  enddo
!$OMP END PARALLEL DO
!      write(*,*)'DM,ibp,xpgrid,ypgrid',20,xpgrid,ypgrid
!      write(*,*)'DM,ibp,FEuler',20,FEuler(xpgrid,ypgrid,:)
!
endsubroutine get_immersed_boundary_rhs
!!-------------------------------------------------
subroutine get_boundary_velocity(ibp,ux,uy,n1,n2)
  implicit none
  integer, intent(in) :: ibp,n1,n2
  double precision,dimension(n1,n2) :: ux,uy
  double precision ::uuxp,uuyp
!
  call from_euler2lag(uuxp,ibp,ux,n1,n2)
  Vxboundary(ibp) = uuxp
  call from_euler2lag(uuyp,ibp,uy,n1,n2)
  Vyboundary(ibp) = uuyp
!  write(*,*) 'DM,ibp,Vxboundary(ibp)',ibp,Vxboundary(ibp),Vyboundary(ibp)
!
end subroutine get_boundary_velocity
!!-------------------------------------------------
subroutine calc_boundary_force_lag(ibp)
  double precision :: xd,yd,displacement
  integer,intent(in) :: ibp
!------
  xd=Xboundary(ibp)-Xboundary_fixed(ibp)
  yd=Yboundary(ibp)-Yboundary_fixed(ibp)
  fboundary(ibp,1)=-elastic_constant*xd*dell(ibp)
  fboundary(ibp,2)=-elastic_constant*yd*dell(ibp)
!------
endsubroutine calc_boundary_force_lag
!!--------------------------------------------------
subroutine from_euler2lag(psip,ibp,psi,n1,n2)
  integer, intent(in) :: n1,n2,ibp
  double precision, dimension(n1,n2), intent(in) :: psi
  double precision, intent(out) :: psip
  integer :: max_neighbour,ineighbour,euler_ix,euler_iy,euler_i,euler_j
  double precision :: eulerx,eulery,rx,ry,psi_ij
  double precision :: xxb,yyb,ddx,ddy
  integer :: xpgrid,ypgrid,xgrid_next,ygrid_next,igrid,jgrid
!---------------
  psip=0.
  xxb=Xboundary(ibp); yyb=Yboundary(ibp)
  call get_grid_loc(xxb,n1,dx,xpgrid,xgrid_next,ddx)
  call get_grid_loc(yyb,n2,dy,ypgrid,ygrid_next,ddy)
  do ineighbour=1,nneighbour
    igrid=neighbour_coord_list(ineighbour,1)
    jgrid=neighbour_coord_list(ineighbour,2)
    euler_ix=xpgrid+igrid; euler_iy=ypgrid+jgrid
    call foldback_pbc(euler_i,euler_ix,n1)
    call foldback_pbc(euler_j,euler_iy,n2)
    psi_ij=psi(euler_i,euler_j)
    rx=abs(igrid*dx-ddx); ry=abs(jgrid*dy-ddy)
    psip=psip+psi_ij*peskin_delta(rx,dx)*peskin_delta(ry,dy)*dx*dy
!    write(*,*)'DM,xpgrid,ypgrid',ibp,xpgrid,ypgrid
!    write(*,*)'DM,ibp,psip',ibp,psip,psi_ij,rx,ry
  enddo
!------------
endsubroutine from_euler2lag
!!-------------------------------------------------
subroutine from_lag2euler(ibp,force_lag)
  integer, intent(in) :: ibp
  double precision,dimension(2), intent(in) :: force_lag
  double precision, dimension(2) :: fij
  integer :: max_neighbour,ineighbour,euler_ix,euler_iy,euler_i,euler_j
  integer :: igrid,jgrid,xgrid_next,ygrid_next,xpgrid,ypgrid
  double precision :: eulerx,eulery,rx,ry,xxb,yyb,ddx,ddy
!--------
  xxb=Xboundary(ibp); yyb=Yboundary(ibp)
  call get_grid_loc(xxb,n1,dx,xpgrid,xgrid_next,ddx)
  call get_grid_loc(yyb,n2,dy,ypgrid,ygrid_next,ddy)
  do ineighbour=1,nneighbour
    igrid=neighbour_coord_list(ineighbour,1)
    jgrid=neighbour_coord_list(ineighbour,2)
    euler_ix=xpgrid+igrid; euler_iy=ypgrid+jgrid
    call foldback_pbc(euler_i,euler_ix,n1)
    call foldback_pbc(euler_j,euler_iy,n2)
    rx=abs(igrid*dx-ddx); ry=abs(jgrid*dy-ddy)
    fij(:) = force_lag(:)*peskin_delta(rx,dx)*peskin_delta(ry,dy)
    FEuler(euler_i,euler_j,1:2) = FEuler(euler_i,euler_j,1:2)+fij(1:2)
!    if (ibp.eq.20) then
!       write(*,*) 'DM ipb,lag force',ibp,force_lag(:)
!       write(*,*) 'DM ineighbour,eul force',ineighbour,fij(:)
!    endif 
  enddo
!  if (ibp.eq.20) then
!    write(*,*) 'DM ibp,eul force',ibp,FEuler(xpgrid,ypgrid,:)
!  endif 
endsubroutine from_lag2euler
!!----------------------------------------------------------
subroutine add_boundary_force(wXu,i1,i2)
  integer,intent(in) :: i1,i2
  double precision,dimension(2) :: wXu 
  if( (i1.eq.138).and.(i2.eq.64)  ) then
!    write(*,*) 'DM wXu before:',wXu
  endif
  wXu = wXu + FEuler(i1,i2,:)
!  if( (i1.eq.138).and.(i2.eq.64)  ) then
!    write(*,*) 'DM i1,i2,FEuler',138,64,FEUler(i1,i2,:)
!    write(*,*) 'DM wXu after:',wXu
!  endif
endsubroutine add_boundary_force
!!----------------------------------------------------------
subroutine evolve_immersed_boundary_euler(dt)
  double precision,intent(in) :: dt
!!
!
!Euler step for the boundary particles. 
!
!$OMP WORKSHARE
  Xboundary=Xboundary+Vxboundary*dt
  Yboundary=Yboundary+Vyboundary*dt
!$OMP END WORKSHARE
  if(lmoving_boundary) then
    call fatal_error('evolve_immersed_boundary_euler','moving boundary not implemented')
!
! When moving boundaries are implemented, chage the position of Xboundary_fixed and
! Yboundary_fixed here. 
!
  endif
!
endsubroutine evolve_immersed_boundary_euler
!-----------------------------------------------
subroutine evolve_immersed_boundary(dt)
  double precision, intent(in) :: dt
  call evolve_immersed_boundary_euler(dt)
endsubroutine evolve_immersed_boundary
!%------------------------------
!
end module immersed_boundary
